{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "view-in-github"
   },
   "source": [
    "<a href=\"https://colab.research.google.com/github/NeuromatchAcademy/course-content/blob/master/tutorials/W3D2_DynamicNetworks/student/W3D2_Tutorial1.ipynb\" target=\"_parent\"><img src=\"https://colab.research.google.com/assets/colab-badge.svg\" alt=\"Open In Colab\"/></a>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text"
   },
   "source": [
    "# Neuromatch Academy: Week 2, Day 4, Tutorial 1\n",
    "# Neuronal Network Dynamics: Neural Rate Models\n",
    "\n",
    "__Content creators:__ Qinglong Gu, Songtin Li, Arvind Kumar, John Murray, Julijana Gjorgjieva \n",
    "\n",
    "__Content reviewers:__  Maryam Vaziri-Pashkam, Ella Batty, Lorenzo Fontolan, Richard Gao, Spiros Chavlis, Michael Waskom\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text"
   },
   "source": [
    "---\n",
    "# Tutorial Objectives\n",
    "\n",
    "The brain is a complex system, not because it is composed of a large number of diverse types of neurons, but mainly because of how neurons are connected to each other. The brain is indeed a network of highly specialized neuronal networks. \n",
    "\n",
    "The activity of a neural network constantly evolves in time. For this reason, neurons can be modeled as dynamical systems. The dynamical system approach is only one of the many modeling approaches that computational neuroscientists have developed (other points of view include information processing,  statistical models, etc.). \n",
    "\n",
    "How the dynamics of neuronal networks affect the representation and processing of information in the brain is an open question. However, signatures of altered brain dynamics present in many brain diseases (e.g., in epilepsy or Parkinson's disease) tell us that it is crucial to study network activity dynamics if we want to understand the brain.\n",
    "\n",
    "In this tutorial, we will simulate and study one of the simplest models of biological neuronal networks. Instead of modeling and simulating individual excitatory neurons (e.g., LIF models that you implemented yesterday), we will treat them as a single homogeneous population and approximate their dynamics using a single one-dimensional equation describing the evolution of their average spiking rate in time.\n",
    "\n",
    "In this tutorial, we will learn how to build a firing rate model of a single population of excitatory neurons. \n",
    "\n",
    "**Steps:**\n",
    "- Write the equation for the firing rate dynamics of a 1D excitatory population.\n",
    "- Visualize the response of the population as a function of parameters such as threshold level and gain, using the frequency-current (F-I) curve.\n",
    "- Numerically simulate the dynamics of the excitatory population and find the fixed points of the system. \n",
    "- Investigate the stability of the fixed points by linearizing the dynamics around them.\n",
    " \n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text"
   },
   "source": [
    "---\n",
    "# Setup"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "cellView": "both",
    "colab": {},
    "colab_type": "code",
    "execution": {
     "iopub.execute_input": "2021-05-25T01:14:09.698750Z",
     "iopub.status.busy": "2021-05-25T01:14:09.698125Z",
     "iopub.status.idle": "2021-05-25T01:14:10.083065Z",
     "shell.execute_reply": "2021-05-25T01:14:10.083635Z"
    }
   },
   "outputs": [],
   "source": [
    "# Imports\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "import scipy.optimize as opt       # root-finding algorithm"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "cellView": "form",
    "colab": {},
    "colab_type": "code",
    "execution": {
     "iopub.execute_input": "2021-05-25T01:14:10.088122Z",
     "iopub.status.busy": "2021-05-25T01:14:10.087529Z",
     "iopub.status.idle": "2021-05-25T01:14:10.169650Z",
     "shell.execute_reply": "2021-05-25T01:14:10.169126Z"
    }
   },
   "outputs": [],
   "source": [
    "# @title Figure Settings\n",
    "import ipywidgets as widgets       # interactive display\n",
    "%config InlineBackend.figure_format = 'retina'\n",
    "plt.style.use(\"https://raw.githubusercontent.com/NeuromatchAcademy/course-content/master/nma.mplstyle\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "cellView": "form",
    "colab": {},
    "colab_type": "code",
    "execution": {
     "iopub.execute_input": "2021-05-25T01:14:10.177126Z",
     "iopub.status.busy": "2021-05-25T01:14:10.176562Z",
     "iopub.status.idle": "2021-05-25T01:14:10.180212Z",
     "shell.execute_reply": "2021-05-25T01:14:10.179757Z"
    }
   },
   "outputs": [],
   "source": [
    "# @title Helper functions\n",
    "\n",
    "\n",
    "def plot_fI(x, f):\n",
    "  plt.figure(figsize=(6, 4))  # plot the figure\n",
    "  plt.plot(x, f, 'k')\n",
    "  plt.xlabel('x (a.u.)', fontsize=14)\n",
    "  plt.ylabel('F(x)', fontsize=14)\n",
    "  plt.show()\n",
    "\n",
    "\n",
    "def plot_dr_r(r, drdt, x_fps=None):\n",
    "  plt.figure()\n",
    "  plt.plot(r, drdt, 'k')\n",
    "  plt.plot(r, 0. * r, 'k--')\n",
    "  if x_fps is not None:\n",
    "    plt.plot(x_fps, np.zeros_like(x_fps), \"ko\", ms=12)\n",
    "  plt.xlabel(r'$r$')\n",
    "  plt.ylabel(r'$\\frac{dr}{dt}$', fontsize=20)\n",
    "  plt.ylim(-0.1, 0.1)\n",
    "\n",
    "\n",
    "def plot_dFdt(x, dFdt):\n",
    "  plt.figure()\n",
    "  plt.plot(x, dFdt, 'r')\n",
    "  plt.xlabel('x (a.u.)', fontsize=14)\n",
    "  plt.ylabel('dF(x)', fontsize=14)\n",
    "  plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text"
   },
   "source": [
    "---\n",
    "# Section 1: Neuronal network dynamics"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "cellView": "form",
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 517
    },
    "colab_type": "code",
    "execution": {
     "iopub.execute_input": "2021-05-25T01:14:10.187743Z",
     "iopub.status.busy": "2021-05-25T01:14:10.187202Z",
     "iopub.status.idle": "2021-05-25T01:14:10.232105Z",
     "shell.execute_reply": "2021-05-25T01:14:10.232811Z"
    },
    "outputId": "83f11632-438e-4703-8853-a526b2e2fd35"
   },
   "outputs": [],
   "source": [
    "# @title Video 1: Dynamic networks\n",
    "from IPython.display import YouTubeVideo\n",
    "video = YouTubeVideo(id=\"p848349hPyw\", width=854, height=480, fs=1)\n",
    "print(\"Video available at https://youtube.com/watch?v=\" + video.id)\n",
    "video"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text"
   },
   "source": [
    "## Section 1.1: Dynamics of a single excitatory population\n",
    "\n",
    "Individual neurons respond by spiking. When we average the spikes of neurons in a population, we can define the average firing activity of the population. In this model, we are interested in how the population-averaged firing varies as a function of time and network parameters. Mathematically, we can describe the firing rate dynamic as:\n",
    "\n",
    "\\begin{align}\n",
    "\\tau \\frac{dr}{dt} &= -r + F(w\\cdot r + I_{\\text{ext}})  \\quad\\qquad (1)\n",
    "\\end{align}\n",
    "\n",
    "$r(t)$ represents the average firing rate of the excitatory population at time $t$, $\\tau$ controls the timescale of the evolution of the average firing rate, $w$ denotes the strength (synaptic weight) of the recurrent input to the population, $I_{\\text{ext}}$ represents the external input, and the transfer function $F(\\cdot)$ (which can be related to f-I curve of individual neurons described in the next sections) represents the population activation function in response to all received inputs.\n",
    "\n",
    "To start building the model, please execute the cell below to initialize the simulation parameters."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "cellView": "form",
    "colab": {},
    "colab_type": "code",
    "execution": {
     "iopub.execute_input": "2021-05-25T01:14:10.238901Z",
     "iopub.status.busy": "2021-05-25T01:14:10.238341Z",
     "iopub.status.idle": "2021-05-25T01:14:10.240359Z",
     "shell.execute_reply": "2021-05-25T01:14:10.240703Z"
    }
   },
   "outputs": [],
   "source": [
    "# @markdown *Execute this cell to set default parameters for a single excitatory population model*\n",
    "\n",
    "\n",
    "def default_pars_single(**kwargs):\n",
    "  pars = {}\n",
    "\n",
    "  # Excitatory parameters\n",
    "  pars['tau'] = 1.     # Timescale of the E population [ms]\n",
    "  pars['a'] = 1.2      # Gain of the E population\n",
    "  pars['theta'] = 2.8  # Threshold of the E population\n",
    "\n",
    "  # Connection strength\n",
    "  pars['w'] = 0.  # E to E, we first set it to 0\n",
    "\n",
    "  # External input\n",
    "  pars['I_ext'] = 0.\n",
    "\n",
    "  # simulation parameters\n",
    "  pars['T'] = 20.       # Total duration of simulation [ms]\n",
    "  pars['dt'] = .1       # Simulation time step [ms]\n",
    "  pars['r_init'] = 0.2  # Initial value of E\n",
    "\n",
    "  # External parameters if any\n",
    "  pars.update(kwargs)\n",
    "\n",
    "  # Vector of discretized time points [ms]\n",
    "  pars['range_t'] = np.arange(0, pars['T'], pars['dt'])\n",
    "\n",
    "  return pars"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text"
   },
   "source": [
    "You can now use:\n",
    "- `pars = default_pars_single()` to get all the parameters, and then you can execute `print(pars)` to check these parameters. \n",
    "- `pars = default_pars_single(T=T_sim, dt=time_step)` to set new simulation time and time step\n",
    "- To update an existing parameter dictionary, use `pars['New_para'] = value`\n",
    "\n",
    "Because `pars` is a dictionary, it can be passed to a function that requires individual parameters as arguments using `my_func(**pars)` syntax."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text"
   },
   "source": [
    "## Section 1.2: F-I curves\n",
    "In electrophysiology, a neuron is often characterized by its spike rate output in response to input currents. This is often called the **F-I** curve, denoting the output spike frequency (**F**) in response to different injected currents (**I**). We estimated this for an LIF neuron in yesterday's tutorial.\n",
    "\n",
    "The transfer function $F(\\cdot)$ in Equation $1$ represents the gain of the population as a function of the total input. The gain is often modeled as a sigmoidal function, i.e., more input drive leads to a nonlinear increase in the population firing rate. The output firing rate will eventually saturate for high input values. \n",
    "\n",
    "A sigmoidal $F(\\cdot)$ is parameterized by its gain $a$ and threshold $\\theta$.\n",
    "\n",
    "$$ F(x;a,\\theta) = \\frac{1}{1+\\text{e}^{-a(x-\\theta)}} - \\frac{1}{1+\\text{e}^{a\\theta}}  \\quad(2)$$\n",
    "\n",
    "The argument $x$ represents the input to the population. Note that the second term is chosen so that $F(0;a,\\theta)=0$.\n",
    "\n",
    "Many other transfer functions (generally monotonic) can be also used. Examples are the rectified linear function $ReLU(x)$ or the hyperbolic tangent $tanh(x)$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text"
   },
   "source": [
    "### Exercise 1: Implement F-I curve \n",
    "\n",
    "Let's first investigate the activation functions before simulating the dynamics of the entire population. \n",
    "\n",
    "In this exercise, you will implement a sigmoidal **F-I** curve or transfer function $F(x)$, with gain $a$ and threshold level $\\theta$ as parameters."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "execution": {
     "iopub.execute_input": "2021-05-25T01:14:10.247603Z",
     "iopub.status.busy": "2021-05-25T01:14:10.246564Z",
     "iopub.status.idle": "2021-05-25T01:14:10.248646Z",
     "shell.execute_reply": "2021-05-25T01:14:10.248157Z"
    }
   },
   "outputs": [],
   "source": [
    "def F(x, a, theta):\n",
    "  \"\"\"\n",
    "  Population activation function.\n",
    "\n",
    "  Args:\n",
    "    x (float): the population input\n",
    "    a (float): the gain of the function\n",
    "    theta (float): the threshold of the function\n",
    "\n",
    "  Returns:\n",
    "    float: the population activation response F(x) for input x\n",
    "  \"\"\"\n",
    "  #################################################\n",
    "  ## TODO for students: compute f = F(x) ##\n",
    "  # Fill out function and remove\n",
    "  raise NotImplementedError(\"Student excercise: implement the f-I function\")\n",
    "  #################################################\n",
    "\n",
    "  # Define the sigmoidal transfer function f = F(x)\n",
    "  f = ...\n",
    "\n",
    "  return f\n",
    "\n",
    "\n",
    "pars = default_pars_single()  # get default parameters\n",
    "x = np.arange(0, 10, .1)      # set the range of input\n",
    "\n",
    "# Uncomment below to test your function\n",
    "# f = F(x, pars['a'], pars['theta'])\n",
    "# plot_fI(x, f)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "cellView": "both",
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 288
    },
    "colab_type": "text",
    "execution": {
     "iopub.execute_input": "2021-05-25T01:14:10.254529Z",
     "iopub.status.busy": "2021-05-25T01:14:10.253973Z",
     "iopub.status.idle": "2021-05-25T01:14:10.439269Z",
     "shell.execute_reply": "2021-05-25T01:14:10.439855Z"
    },
    "outputId": "23aa765f-8d7d-481d-a204-985e5b7ec535"
   },
   "source": [
    "[*Click for solution*](https://github.com/NeuromatchAcademy/course-content/tree/master//tutorials/W2D4_DynamicNetworks/solutions/W2D4_Tutorial1_Solution_45ddc05f.py)\n",
    "\n",
    "*Example output:*\n",
    "\n",
    "<img alt='Solution hint' align='left' width=416 height=272 src=https://raw.githubusercontent.com/NeuromatchAcademy/course-content/master/tutorials/W2D4_DynamicNetworks/static/W2D4_Tutorial1_Solution_45ddc05f_0.png>\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text"
   },
   "source": [
    "### Interactive Demo: Parameter exploration of F-I curve\n",
    "Here's an interactive demo that shows how the F-I curve changes for different values of the gain and threshold parameters. How do the gain and threshold parameters affect the F-I curve?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "cellView": "form",
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 494,
     "referenced_widgets": [
      "b61bd80d7d364ad7bf04d01c1cd74910",
      "245b4a9a358641198e846046ca13cefa",
      "144fe28586da4c43a999d02b60c7df79",
      "03250c609c324c2c891bf4a4ee465a26",
      "85a36101c38c4aa48deef2ab8946d246",
      "67d602777c7946e0abeb5955b326ce60",
      "95b37093de5d46c981097831e3419405",
      "ee845cd049f94bb19ae8222865495d38",
      "271c3bd6914c4808aa0aa39d6239a80b",
      "a53e1476f34a4e1aae19c69a62e60dc1"
     ]
    },
    "colab_type": "code",
    "execution": {
     "iopub.execute_input": "2021-05-25T01:14:10.469443Z",
     "iopub.status.busy": "2021-05-25T01:14:10.468641Z",
     "iopub.status.idle": "2021-05-25T01:14:10.663111Z",
     "shell.execute_reply": "2021-05-25T01:14:10.663568Z"
    },
    "outputId": "47e4f04f-041d-4573-8114-ba644cb8e313"
   },
   "outputs": [],
   "source": [
    "# @title\n",
    "\n",
    "# @markdown Make sure you execute this cell to enable the widget!\n",
    "\n",
    "\n",
    "def interactive_plot_FI(a, theta):\n",
    "  \"\"\"\n",
    "  Population activation function.\n",
    "\n",
    "  Expecxts:\n",
    "    a     : the gain of the function\n",
    "    theta : the threshold of the function\n",
    "\n",
    "  Returns:\n",
    "    plot the F-I curve with give parameters\n",
    "  \"\"\"\n",
    "\n",
    "  # set the range of input\n",
    "  x = np.arange(0, 10, .1)\n",
    "  plt.figure()\n",
    "  plt.plot(x, F(x, a, theta), 'k')\n",
    "  plt.xlabel('x (a.u.)', fontsize=14)\n",
    "  plt.ylabel('F(x)', fontsize=14)\n",
    "  plt.show()\n",
    "\n",
    "\n",
    "_ = widgets.interact(interactive_plot_FI, a=(0.3, 3, 0.3), theta=(2, 4, 0.2))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab": {},
    "colab_type": "text",
    "execution": {
     "iopub.execute_input": "2021-05-25T01:14:10.668988Z",
     "iopub.status.busy": "2021-05-25T01:14:10.668438Z",
     "iopub.status.idle": "2021-05-25T01:14:10.673261Z",
     "shell.execute_reply": "2021-05-25T01:14:10.672739Z"
    }
   },
   "source": [
    "[*Click for solution*](https://github.com/NeuromatchAcademy/course-content/tree/master//tutorials/W2D4_DynamicNetworks/solutions/W2D4_Tutorial1_Solution_1c0165d7.py)\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text"
   },
   "source": [
    "## Section 1.3: Simulation scheme of E dynamics\n",
    "\n",
    "Because $F(\\cdot)$ is a nonlinear function, the exact solution of Equation $1$ can not be determined via analytical methods. Therefore, numerical methods must be used to find the solution. In practice, the derivative on the left-hand side of Equation $1$ can be approximated using the Euler method on a time-grid of stepsize $\\Delta t$:\n",
    "\n",
    "\\begin{align}\n",
    "&\\frac{dr}{dt} \\approx \\frac{r[k+1]-r[k]}{\\Delta t} \n",
    "\\end{align}\n",
    "where $r[k] = r(k\\Delta t)$. \n",
    "\n",
    "Thus,\n",
    "\n",
    "$$\\Delta r[k] = \\frac{\\Delta t}{\\tau}[-r[k] + F(w\\cdot r[k] + I_{\\text{ext}}[k];a,\\theta)]$$\n",
    "\n",
    "\n",
    "Hence, Equation (1) is updated at each time step by:\n",
    "\n",
    "$$r[k+1] = r[k] + \\Delta r[k]$$\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "cellView": "form",
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 283
    },
    "colab_type": "code",
    "execution": {
     "iopub.execute_input": "2021-05-25T01:14:10.680469Z",
     "iopub.status.busy": "2021-05-25T01:14:10.679879Z",
     "iopub.status.idle": "2021-05-25T01:14:10.684184Z",
     "shell.execute_reply": "2021-05-25T01:14:10.683699Z"
    },
    "outputId": "39ac5bc8-5083-44c1-be6e-2d91c80e9983"
   },
   "outputs": [],
   "source": [
    "# @markdown *Execute this cell to enable the single population rate model simulator: `simulate_single`*\n",
    "\n",
    "\n",
    "def simulate_single(pars):\n",
    "  \"\"\"\n",
    "  Simulate an excitatory population of neurons\n",
    "\n",
    "  Args:\n",
    "    pars   : Parameter dictionary\n",
    "\n",
    "  Returns:\n",
    "    rE     : Activity of excitatory population (array)\n",
    "\n",
    "  Example:\n",
    "    pars = default_pars_single()\n",
    "    r = simulate_single(pars)\n",
    "  \"\"\"\n",
    "\n",
    "  # Set parameters\n",
    "  tau, a, theta = pars['tau'], pars['a'], pars['theta']\n",
    "  w = pars['w']\n",
    "  I_ext = pars['I_ext']\n",
    "  r_init = pars['r_init']\n",
    "  dt, range_t = pars['dt'], pars['range_t']\n",
    "  Lt = range_t.size\n",
    "\n",
    "  # Initialize activity\n",
    "  r = np.zeros(Lt)\n",
    "  r[0] = r_init\n",
    "  I_ext = I_ext * np.ones(Lt)\n",
    "\n",
    "  # Update the E activity\n",
    "  for k in range(Lt - 1):\n",
    "      dr = dt / tau * (-r[k] + F(w * r[k] + I_ext[k], a, theta))\n",
    "      r[k+1] = r[k] + dr\n",
    "\n",
    "  return r\n",
    "\n",
    "help(simulate_single)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text"
   },
   "source": [
    "### Interactive Demo: Parameter Exploration of single population dynamics\n",
    "\n",
    "Note that $w=0$, as in the default setting, means no recurrent input to the neuron population in Equation (1). Hence, the dynamics are entirely determined by the external input $I_{\\text{ext}}$. Explore these dynamics in this interactive demo.\n",
    "\n",
    "How does $r_{\\text{sim}}(t)$ change with different $I_{\\text{ext}}$ values? How does it change with different $\\tau$ values? Investigate the relationship between $F(I_{\\text{ext}}; a, \\theta)$ and the steady value of $r(t)$. \n",
    "\n",
    "Note that, $r_{\\rm ana}(t)$ denotes the analytical solution - you will learn how this is computed in the next section."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "cellView": "form",
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 494,
     "referenced_widgets": [
      "4701ae22a90d47b49c0102c1d74dea93",
      "724ab274c39d4e00ac6ea1a63a3e5d29",
      "d3012a33cb774402874129a452dbc74e",
      "f276d852e7cc4313a129277277081c91",
      "80f7139b2cc04498ab0157bb9bfe400f",
      "808a66d0557e4994aa6098643e9b1994",
      "fd5bc6df0676470bb62e6dadf3ef1aae",
      "c62ce0f30fc34364a816de5affce9108",
      "79539b376eb1454dbcec5b8a0bc058b1",
      "85229a0353be4da285a280bff960af5f"
     ]
    },
    "colab_type": "code",
    "execution": {
     "iopub.execute_input": "2021-05-25T01:14:10.715350Z",
     "iopub.status.busy": "2021-05-25T01:14:10.695091Z",
     "iopub.status.idle": "2021-05-25T01:14:11.192495Z",
     "shell.execute_reply": "2021-05-25T01:14:11.191951Z"
    },
    "outputId": "91dd1a42-b4ed-400d-e6e3-7a5186fec0b4"
   },
   "outputs": [],
   "source": [
    "# @title\n",
    "\n",
    "# @markdown Make sure you execute this cell to enable the widget!\n",
    "\n",
    "# get default parameters\n",
    "pars = default_pars_single(T=20.)\n",
    "\n",
    "\n",
    "def Myplot_E_diffI_difftau(I_ext, tau):\n",
    "  # set external input and time constant\n",
    "  pars['I_ext'] = I_ext\n",
    "  pars['tau'] = tau\n",
    "\n",
    "  # simulation\n",
    "  r = simulate_single(pars)\n",
    "\n",
    "  # Analytical Solution\n",
    "  r_ana = (pars['r_init']\n",
    "           + (F(I_ext, pars['a'], pars['theta'])\n",
    "           - pars['r_init']) * (1. - np.exp(-pars['range_t'] / pars['tau'])))\n",
    "\n",
    "  # plot\n",
    "  plt.figure()\n",
    "  plt.plot(pars['range_t'], r, 'b', label=r'$r_{\\mathrm{sim}}$(t)', alpha=0.5,\n",
    "           zorder=1)\n",
    "  plt.plot(pars['range_t'], r_ana, 'b--', lw=5, dashes=(2, 2),\n",
    "           label=r'$r_{\\mathrm{ana}}$(t)', zorder=2)\n",
    "  plt.plot(pars['range_t'],\n",
    "           F(I_ext, pars['a'], pars['theta']) * np.ones(pars['range_t'].size),\n",
    "           'k--', label=r'$F(I_{\\mathrm{ext}})$')\n",
    "  plt.xlabel('t (ms)', fontsize=16.)\n",
    "  plt.ylabel('Activity r(t)', fontsize=16.)\n",
    "  plt.legend(loc='best', fontsize=14.)\n",
    "  plt.show()\n",
    "\n",
    "\n",
    "_ = widgets.interact(Myplot_E_diffI_difftau, I_ext=(0.0, 10., 1.),\n",
    "                     tau=(1., 5., 0.2))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab": {},
    "colab_type": "text",
    "execution": {
     "iopub.execute_input": "2021-05-25T01:14:11.199139Z",
     "iopub.status.busy": "2021-05-25T01:14:11.198599Z",
     "iopub.status.idle": "2021-05-25T01:14:11.201275Z",
     "shell.execute_reply": "2021-05-25T01:14:11.201727Z"
    }
   },
   "source": [
    "[*Click for solution*](https://github.com/NeuromatchAcademy/course-content/tree/master//tutorials/W2D4_DynamicNetworks/solutions/W2D4_Tutorial1_Solution_65dee3e7.py)\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text"
   },
   "source": [
    "## Think!\n",
    "Above, we have numerically solved a system driven by a positive input. Yet, $r_E(t)$ either decays to zero or reaches a fixed non-zero value.\n",
    "- Why doesn't the solution of the system \"explode\" in a finite time? In other words, what guarantees that $r_E$(t) stays finite? \n",
    "- Which parameter would you change in order to increase the maximum value of the response? "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab": {},
    "colab_type": "text",
    "execution": {
     "iopub.execute_input": "2021-05-25T01:14:11.206313Z",
     "iopub.status.busy": "2021-05-25T01:14:11.205745Z",
     "iopub.status.idle": "2021-05-25T01:14:11.208396Z",
     "shell.execute_reply": "2021-05-25T01:14:11.208816Z"
    }
   },
   "source": [
    "[*Click for solution*](https://github.com/NeuromatchAcademy/course-content/tree/master//tutorials/W2D4_DynamicNetworks/solutions/W2D4_Tutorial1_Solution_5a95a98e.py)\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text"
   },
   "source": [
    "---\n",
    "# Section 2: Fixed points of the single population system\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "cellView": "form",
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 517
    },
    "colab_type": "code",
    "execution": {
     "iopub.execute_input": "2021-05-25T01:14:11.214334Z",
     "iopub.status.busy": "2021-05-25T01:14:11.213808Z",
     "iopub.status.idle": "2021-05-25T01:14:11.257011Z",
     "shell.execute_reply": "2021-05-25T01:14:11.257435Z"
    },
    "outputId": "f25d4b12-6499-414e-bee7-88ea1ea9c557"
   },
   "outputs": [],
   "source": [
    "# @title Video 2: Fixed point\n",
    "from IPython.display import YouTubeVideo\n",
    "video = YouTubeVideo(id=\"Ox3ELd1UFyo\", width=854, height=480, fs=1)\n",
    "print(\"Video available at https://youtube.com/watch?v=\" + video.id)\n",
    "video"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text"
   },
   "source": [
    "As you varied the two parameters in the last Interactive Demo, you noticed that, while at first the system output quickly changes, with time, it reaches its maximum/minimum value and does not change anymore. The value eventually reached by the system is called the **steady state** of the system, or the **fixed point**. Essentially, in the steady states the derivative with respect to time of the activity ($r$) is zero, i.e. $\\displaystyle \\frac{dr}{dt}=0$. \n",
    "\n",
    "We can find that the steady state of the Equation. (1) by setting $\\displaystyle{\\frac{dr}{dt}=0}$ and solve for $r$:\n",
    "\n",
    "$$-r_{\\text{steady}} + F(w\\cdot r_{\\text{steady}} + I_{\\text{ext}};a,\\theta) = 0, \\qquad (3)$$\n",
    "\n",
    "When it exists, the solution of Equation. (3) defines a **fixed point** of the dynamical system in Equation (1). Note that if $F(x)$ is nonlinear, it is not always possible to find an analytical solution, but the solution can be found via numerical simulations, as we will do later.\n",
    "\n",
    "From the Interactive Demo, one could also notice that the value of $\\tau$ influences how quickly the activity will converge to the steady state from its initial value. \n",
    "\n",
    "In the specific case of $w=0$, we can also analytically compute  the solution of Equation (1) (i.e., the thick blue dashed line) and deduce the role of $\\tau$ in determining the convergence to the fixed point: \n",
    "\n",
    "$$\\displaystyle{r(t) = \\big{[}F(I_{\\text{ext}};a,\\theta) -r(t=0)\\big{]} (1-\\text{e}^{-\\frac{t}{\\tau}})} + r(t=0)$$ \\\\\n",
    "\n",
    "We can now numerically calculate the fixed point with a root finding algorithm."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text"
   },
   "source": [
    "## Exercise 2: Visualization of the fixed points\n",
    "\n",
    "When it is not possible to find the solution for Equation (3) analytically, a graphical approach can be taken. To that end, it is useful to plot $\\displaystyle{\\frac{dr}{dt}}$ as a function of $r$. The values of $r$ for which the plotted function crosses zero on the y axis correspond to fixed points. \n",
    "\n",
    "Here, let us, for example, set $w=5.0$ and $I^{\\text{ext}}=0.5$. From Equation (1), you can obtain\n",
    "\n",
    "$$\\frac{dr}{dt} = [-r + F(w\\cdot r + I^{\\text{ext}})]\\,/\\,\\tau $$\n",
    "\n",
    "Then, plot the $dr/dt$ as a function of $r$, and check for the presence of fixed points. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "execution": {
     "iopub.execute_input": "2021-05-25T01:14:11.263163Z",
     "iopub.status.busy": "2021-05-25T01:14:11.261994Z",
     "iopub.status.idle": "2021-05-25T01:14:11.263779Z",
     "shell.execute_reply": "2021-05-25T01:14:11.264202Z"
    }
   },
   "outputs": [],
   "source": [
    "def compute_drdt(r, I_ext, w, a, theta, tau, **other_pars):\n",
    "  \"\"\"Given parameters, compute dr/dt as a function of r.\n",
    "\n",
    "  Args:\n",
    "    r (1D array) : Average firing rate of the excitatory population\n",
    "    I_ext, w, a, theta, tau (numbers): Simulation parameters to use\n",
    "    other_pars : Other simulation parameters are unused by this function\n",
    "\n",
    "  Returns\n",
    "    drdt function for each value of r\n",
    "  \"\"\"\n",
    "  #########################################################################\n",
    "  # TODO compute drdt and disable the error\n",
    "  raise NotImplementedError(\"Finish the compute_drdt function\")\n",
    "  #########################################################################\n",
    "\n",
    "  # Calculate drdt\n",
    "  drdt = ...\n",
    "\n",
    "  return drdt\n",
    "\n",
    "\n",
    "# Define a vector of r values and the simulation parameters\n",
    "r = np.linspace(0, 1, 1000)\n",
    "pars = default_pars_single(I_ext=0.5, w=5)\n",
    "\n",
    "# Uncomment to test your function\n",
    "# drdt = compute_drdt(r, **pars)\n",
    "# plot_dr_r(r, drdt)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 485
    },
    "colab_type": "text",
    "execution": {
     "iopub.execute_input": "2021-05-25T01:14:11.269164Z",
     "iopub.status.busy": "2021-05-25T01:14:11.268646Z",
     "iopub.status.idle": "2021-05-25T01:14:11.544368Z",
     "shell.execute_reply": "2021-05-25T01:14:11.543854Z"
    },
    "outputId": "406680a9-25af-446f-df62-6b7ba39b4e12"
   },
   "source": [
    "[*Click for solution*](https://github.com/NeuromatchAcademy/course-content/tree/master//tutorials/W2D4_DynamicNetworks/solutions/W2D4_Tutorial1_Solution_c5280901.py)\n",
    "\n",
    "*Example output:*\n",
    "\n",
    "<img alt='Solution hint' align='left' width=558 height=413 src=https://raw.githubusercontent.com/NeuromatchAcademy/course-content/master/tutorials/W2D4_DynamicNetworks/static/W2D4_Tutorial1_Solution_c5280901_0.png>\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text"
   },
   "source": [
    "## Exercise 3: Fixed point calculation\n",
    "\n",
    "We will now find the fixed points numerically. To do so, we need to specif initial values ($r_{\\text{guess}}$) for the root-finding algorithm to start from. From the line $\\displaystyle{\\frac{dr}{dt}}$ plotted above in Exercise 2, initial values can be chosen as a set of values close to where the line crosses zero on the y axis (real fixed point).\n",
    "\n",
    "The next cell defines three helper functions that we will use:\n",
    "\n",
    "- `my_fp_single(r_guess, **pars)` uses a root-finding algorithm to locate a fixed point near a given initial value\n",
    "- `check_fp_single(x_fp, **pars)`, verifies that the values of $r_{\\rm fp}$ for which $\\displaystyle{\\frac{dr}{dt}} = 0$ are the true fixed points\n",
    "- `my_fp_finder(r_guess_vector, **pars)` accepts an array of initial values and finds the same number of fixed points, using the above two functions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "cellView": "form",
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 247
    },
    "colab_type": "code",
    "execution": {
     "iopub.execute_input": "2021-05-25T01:14:11.552991Z",
     "iopub.status.busy": "2021-05-25T01:14:11.552447Z",
     "iopub.status.idle": "2021-05-25T01:14:11.554212Z",
     "shell.execute_reply": "2021-05-25T01:14:11.554630Z"
    },
    "outputId": "079fa0d8-885f-4f21-eca5-6904468e68ee"
   },
   "outputs": [],
   "source": [
    "# @markdown *Execute this cell to enable the fixed point functions*\n",
    "\n",
    "def my_fp_single(r_guess, a, theta, w, I_ext, **other_pars):\n",
    "  \"\"\"\n",
    "  Calculate the fixed point through drE/dt=0\n",
    "\n",
    "  Args:\n",
    "    r_guess  : Initial value used for scipy.optimize function\n",
    "    a, theta, w, I_ext : simulation parameters\n",
    "\n",
    "  Returns:\n",
    "    x_fp    : value of fixed point\n",
    "  \"\"\"\n",
    "  # define the right hand of E dynamics\n",
    "  def my_WCr(x):\n",
    "    r = x\n",
    "    drdt = (-r + F(w * r + I_ext, a, theta))\n",
    "    y = np.array(drdt)\n",
    "\n",
    "    return y\n",
    "\n",
    "  x0 = np.array(r_guess)\n",
    "  x_fp = opt.root(my_WCr, x0).x.item()\n",
    "\n",
    "  return x_fp\n",
    "\n",
    "\n",
    "def check_fp_single(x_fp, a, theta, w, I_ext, mytol=1e-4, **other_pars):\n",
    "  \"\"\"\n",
    "   Verify |dr/dt| < mytol\n",
    "\n",
    "  Args:\n",
    "    fp      : value of fixed point\n",
    "    a, theta, w, I_ext: simulation parameters\n",
    "    mytol   : tolerance, default as 10^{-4}\n",
    "\n",
    "  Returns :\n",
    "    Whether it is a correct fixed point: True/False\n",
    "  \"\"\"\n",
    "  # calculate Equation(3)\n",
    "  y = x_fp - F(w * x_fp + I_ext, a, theta)\n",
    "\n",
    "  # Here we set tolerance as 10^{-4}\n",
    "  return np.abs(y) < mytol\n",
    "\n",
    "\n",
    "def my_fp_finder(pars, r_guess_vector, mytol=1e-4):\n",
    "  \"\"\"\n",
    "  Calculate the fixed point(s) through drE/dt=0\n",
    "\n",
    "  Args:\n",
    "    pars    : Parameter dictionary\n",
    "    r_guess_vector  : Initial values used for scipy.optimize function\n",
    "    mytol   : tolerance for checking fixed point, default as 10^{-4}\n",
    "\n",
    "  Returns:\n",
    "    x_fps   : values of fixed points\n",
    "\n",
    "  \"\"\"\n",
    "  x_fps = []\n",
    "  correct_fps = []\n",
    "  for r_guess in r_guess_vector:\n",
    "    x_fp = my_fp_single(r_guess, **pars)\n",
    "    if check_fp_single(x_fp, **pars, mytol=mytol):\n",
    "      x_fps.append(x_fp)\n",
    "\n",
    "  return x_fps\n",
    "\n",
    "help(my_fp_finder)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "execution": {
     "iopub.execute_input": "2021-05-25T01:14:11.563327Z",
     "iopub.status.busy": "2021-05-25T01:14:11.561813Z",
     "iopub.status.idle": "2021-05-25T01:14:11.565118Z",
     "shell.execute_reply": "2021-05-25T01:14:11.564279Z"
    }
   },
   "outputs": [],
   "source": [
    "r = np.linspace(0, 1, 1000)\n",
    "pars = default_pars_single(I_ext=0.5, w=5)\n",
    "drdt = compute_drdt(r, **pars)\n",
    "\n",
    "#############################################################################\n",
    "# TODO for students:\n",
    "# Define initial values close to the intersections of drdt and y=0\n",
    "# (How many initial values? Hint: How many times do the two lines intersect?)\n",
    "# Calculate the fixed point with these initial values and plot them\n",
    "#############################################################################\n",
    "r_guess_vector = [...]\n",
    "\n",
    "# Uncomment to test your values\n",
    "# x_fps = my_fp_finder(pars, r_guess_vector)\n",
    "# plot_dr_r(r, drdt, x_fps)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 485
    },
    "colab_type": "text",
    "execution": {
     "iopub.execute_input": "2021-05-25T01:14:11.584800Z",
     "iopub.status.busy": "2021-05-25T01:14:11.584190Z",
     "iopub.status.idle": "2021-05-25T01:14:11.853546Z",
     "shell.execute_reply": "2021-05-25T01:14:11.852289Z"
    },
    "outputId": "4cf2cddb-0cc0-45be-e577-48e69a6ee6f3"
   },
   "source": [
    "[*Click for solution*](https://github.com/NeuromatchAcademy/course-content/tree/master//tutorials/W2D4_DynamicNetworks/solutions/W2D4_Tutorial1_Solution_0637b6bf.py)\n",
    "\n",
    "*Example output:*\n",
    "\n",
    "<img alt='Solution hint' align='left' width=558 height=413 src=https://raw.githubusercontent.com/NeuromatchAcademy/course-content/master/tutorials/W2D4_DynamicNetworks/static/W2D4_Tutorial1_Solution_0637b6bf_0.png>\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text"
   },
   "source": [
    "## Interactive Demo: fixed points as a function of recurrent and external inputs.\n",
    "\n",
    "You can now explore how the previous plot changes when the recurrent coupling $w$ and the external input $I_{\\text{ext}}$ take different values. How does the number of fixed points change?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "cellView": "form",
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 494,
     "referenced_widgets": [
      "875329f08b434412a4e553ad777cf762",
      "abb9d829af484d84be1463074b52f5b4",
      "04f4018e9c2544618446287713fc8876",
      "b4066b09119c4faca10ee55f9438f63b",
      "b07bc0d64ad5444cb239d0ac25bc164f",
      "0e1d4030d1c44fdab228e94b40f5a66d",
      "1a161bda1e3e405289c0fc5d1f927aa9",
      "1e1c2f33a69842cb86d1a931e41182ad",
      "12a33e3c45bf43bf9e64366527bba3f6",
      "408c7f906c60459fb4375e995b8064c2"
     ]
    },
    "colab_type": "code",
    "execution": {
     "iopub.execute_input": "2021-05-25T01:14:11.884364Z",
     "iopub.status.busy": "2021-05-25T01:14:11.883747Z",
     "iopub.status.idle": "2021-05-25T01:14:12.124499Z",
     "shell.execute_reply": "2021-05-25T01:14:12.124000Z"
    },
    "outputId": "b87c83c5-d95c-40a2-efa5-e657bbec0f2a"
   },
   "outputs": [],
   "source": [
    "# @title\n",
    "\n",
    "# @markdown Make sure you execute this cell to enable the widget!\n",
    "\n",
    "\n",
    "def plot_intersection_single(w, I_ext):\n",
    "  # set your parameters\n",
    "  pars = default_pars_single(w=w, I_ext=I_ext)\n",
    "\n",
    "  # find fixed points\n",
    "  r_init_vector = [0, .4, .9]\n",
    "  x_fps = my_fp_finder(pars, r_init_vector)\n",
    "\n",
    "  # plot\n",
    "  r = np.linspace(0, 1., 1000)\n",
    "  drdt = (-r + F(w * r + I_ext, pars['a'], pars['theta'])) / pars['tau']\n",
    "\n",
    "  plot_dr_r(r, drdt, x_fps)\n",
    "\n",
    "_ = widgets.interact(plot_intersection_single, w=(1, 7, 0.2),\n",
    "                     I_ext=(0, 3, 0.1))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab": {},
    "colab_type": "text",
    "execution": {
     "iopub.execute_input": "2021-05-25T01:14:12.130652Z",
     "iopub.status.busy": "2021-05-25T01:14:12.130113Z",
     "iopub.status.idle": "2021-05-25T01:14:12.132637Z",
     "shell.execute_reply": "2021-05-25T01:14:12.133056Z"
    }
   },
   "source": [
    "[*Click for solution*](https://github.com/NeuromatchAcademy/course-content/tree/master//tutorials/W2D4_DynamicNetworks/solutions/W2D4_Tutorial1_Solution_20486792.py)\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text"
   },
   "source": [
    "---\n",
    "# Summary\n",
    "\n",
    "In this tutorial, we have investigated the dynamics of a rate-based single population of neurons.\n",
    "\n",
    "We learned about:\n",
    "- The effect of the input parameters and the time constant of the network on the dynamics of the population.\n",
    "- How to find the fixed point(s) of the system.\n",
    "\n",
    "Next, we have two Bonus, but important concepts in dynamical system analysis and simulation. If you have time left, watch the next video and proceed to solve the exercises. You will learn:\n",
    "\n",
    "- How to determine the stability of a fixed point by linearizing the system.\n",
    "- How to add realistic inputs to our model."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text"
   },
   "source": [
    "---\n",
    "# Bonus 1: Stability of a fixed point"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "cellView": "form",
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 517
    },
    "colab_type": "code",
    "execution": {
     "iopub.execute_input": "2021-05-25T01:14:12.138286Z",
     "iopub.status.busy": "2021-05-25T01:14:12.137713Z",
     "iopub.status.idle": "2021-05-25T01:14:12.180907Z",
     "shell.execute_reply": "2021-05-25T01:14:12.181334Z"
    },
    "outputId": "11b2ff78-d268-404b-8b48-596a70315854"
   },
   "outputs": [],
   "source": [
    "# @title Video 3: Stability of fixed points\n",
    "from IPython.display import YouTubeVideo\n",
    "video = YouTubeVideo(id=\"KKMlWWU83Jg\", width=854, height=480, fs=1)\n",
    "print(\"Video available at https://youtube.com/watch?v=\" + video.id)\n",
    "video"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text"
   },
   "source": [
    "#### Initial values and trajectories\n",
    "\n",
    "Here, let us first set $w=5.0$ and $I_{\\text{ext}}=0.5$, and investigate the dynamics of $r(t)$ starting with different initial values $r(0) \\equiv r_{\\text{init}}$. We will plot the trajectories of $r(t)$ with $r_{\\text{init}} = 0.0, 0.1, 0.2,..., 0.9$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "cellView": "form",
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 357
    },
    "colab_type": "code",
    "execution": {
     "iopub.execute_input": "2021-05-25T01:14:12.203618Z",
     "iopub.status.busy": "2021-05-25T01:14:12.201252Z",
     "iopub.status.idle": "2021-05-25T01:14:12.658831Z",
     "shell.execute_reply": "2021-05-25T01:14:12.659300Z"
    },
    "outputId": "e69ca3c9-9981-4fd4-e9bf-73d62cf764de"
   },
   "outputs": [],
   "source": [
    "# @markdown Execute this cell to see the trajectories!\n",
    "\n",
    "pars = default_pars_single()\n",
    "pars['w'] = 5.0\n",
    "pars['I_ext'] = 0.5\n",
    "\n",
    "plt.figure(figsize=(8, 5))\n",
    "for ie in range(10):\n",
    "  pars['r_init'] = 0.1 * ie  # set the initial value\n",
    "  r = simulate_single(pars)  # run the simulation\n",
    "\n",
    "  # plot the activity with given initial\n",
    "  plt.plot(pars['range_t'], r, 'b', alpha=0.1 + 0.1 * ie,\n",
    "           label=r'r$_{\\mathrm{init}}$=%.1f' % (0.1 * ie))\n",
    "\n",
    "plt.xlabel('t (ms)')\n",
    "plt.title('Two steady states?')\n",
    "plt.ylabel(r'$r$(t)')\n",
    "plt.legend(loc=[1.01, -0.06], fontsize=14)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text"
   },
   "source": [
    "## Interactive Demo: dynamics as a function of the initial value\n",
    "\n",
    "Let's now set $r_{\\rm init}$ to a value of your choice in this demo. How does the solution change? What do you observe?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "cellView": "form",
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 462,
     "referenced_widgets": [
      "a861ef888bd041ad8c47551e665cc9ba",
      "8d37cd431d964d96a50582482ff028d4",
      "ff23d44495b44d729861ba6cf52c0bca",
      "ba13f61549f24b14ba70e387aa1a9485",
      "f9e0df8ceaf647949c1c8975651dfda2",
      "cbec43cfc0064abf805c1d52bb395caa",
      "ed7948064a3b415d8429da73484125f7"
     ]
    },
    "colab_type": "code",
    "execution": {
     "iopub.execute_input": "2021-05-25T01:14:12.683404Z",
     "iopub.status.busy": "2021-05-25T01:14:12.682769Z",
     "iopub.status.idle": "2021-05-25T01:14:12.872627Z",
     "shell.execute_reply": "2021-05-25T01:14:12.873079Z"
    },
    "outputId": "2962d6fb-cad4-45d2-8b39-910437af7b8b"
   },
   "outputs": [],
   "source": [
    "# @title\n",
    "\n",
    "# @markdown Make sure you execute this cell to enable the widget!\n",
    "\n",
    "pars = default_pars_single(w=5.0, I_ext=0.5)\n",
    "\n",
    "def plot_single_diffEinit(r_init):\n",
    "  pars['r_init'] = r_init\n",
    "  r = simulate_single(pars)\n",
    "\n",
    "  plt.figure()\n",
    "  plt.plot(pars['range_t'], r, 'b', zorder=1)\n",
    "  plt.plot(0, r[0], 'bo', alpha=0.7, zorder=2)\n",
    "  plt.xlabel('t (ms)', fontsize=16)\n",
    "  plt.ylabel(r'$r(t)$', fontsize=16)\n",
    "  plt.ylim(0, 1.0)\n",
    "  plt.show()\n",
    "\n",
    "\n",
    "_ = widgets.interact(plot_single_diffEinit, r_init=(0, 1, 0.02))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab": {},
    "colab_type": "text",
    "execution": {
     "iopub.execute_input": "2021-05-25T01:14:12.875911Z",
     "iopub.status.busy": "2021-05-25T01:14:12.875392Z",
     "iopub.status.idle": "2021-05-25T01:14:12.879393Z",
     "shell.execute_reply": "2021-05-25T01:14:12.878963Z"
    }
   },
   "source": [
    "[*Click for solution*](https://github.com/NeuromatchAcademy/course-content/tree/master//tutorials/W2D4_DynamicNetworks/solutions/W2D4_Tutorial1_Solution_4d2de6a0.py)\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text"
   },
   "source": [
    "### Stability analysis via linearization of the dynamics\n",
    "\n",
    "Just like Equation $1$ in the case ($w=0$) discussed above, a generic linear system \n",
    "$$\\frac{dx}{dt} = \\lambda (x - b),$$ \n",
    "has a fixed point for $x=b$. The analytical solution of such a system can be found to be:\n",
    "$$x(t) = b + \\big{(} x(0) - b \\big{)} \\text{e}^{\\lambda t}.$$ \n",
    "Now consider a small perturbation of the activity around the fixed point: $x(0) = b+ \\epsilon$, where $|\\epsilon| \\ll 1$. Will the perturbation $\\epsilon(t)$ grow with time or will it decay to the fixed point? The evolution of the perturbation with time can be written, using the analytical solution for $x(t)$, as:\n",
    " $$\\epsilon (t) = x(t) - b = \\epsilon \\text{e}^{\\lambda t}$$\n",
    "\n",
    "- if $\\lambda < 0$, $\\epsilon(t)$ decays to zero, $x(t)$ will still converge to $b$ and the fixed point is \"**stable**\".\n",
    "\n",
    "- if $\\lambda > 0$, $\\epsilon(t)$ grows with time, $x(t)$ will leave the fixed point $b$ exponentially, and the fixed point is, therefore, \"**unstable**\" ."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text"
   },
   "source": [
    "### Compute the stability of Equation $1$\n",
    "\n",
    "Similar to what we did in the linear system above, in order to determine the stability of a fixed point $r^{*}$ of the excitatory population dynamics, we perturb Equation (1) around $r^{*}$ by $\\epsilon$, i.e. $r = r^{*} + \\epsilon$. We can plug in Equation (1) and obtain the equation determining the time evolution of the perturbation $\\epsilon(t)$:\n",
    "\n",
    "\\begin{align}\n",
    "\\tau \\frac{d\\epsilon}{dt} \\approx -\\epsilon + w F'(w\\cdot r^{*} + I_{\\text{ext}};a,\\theta) \\epsilon \n",
    "\\end{align}\n",
    "\n",
    "where $F'(\\cdot)$ is the derivative of the transfer function $F(\\cdot)$. We can rewrite the above equation as:\n",
    "\n",
    "\\begin{align}\n",
    "\\frac{d\\epsilon}{dt} \\approx \\frac{\\epsilon}{\\tau }[-1 + w F'(w\\cdot r^* + I_{\\text{ext}};a,\\theta)] \n",
    "\\end{align}\n",
    "\n",
    "That is, as in the linear system above, the value of\n",
    "\n",
    "$$\\lambda = [-1+ wF'(w\\cdot r^* + I_{\\text{ext}};a,\\theta)]/\\tau \\qquad (4)$$\n",
    "\n",
    "determines whether the perturbation will grow or decay to zero, i.e., $\\lambda$ defines the stability of the fixed point. This value is called the **eigenvalue** of the dynamical system."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text"
   },
   "source": [
    "## Exercise 4: Compute $dF$\n",
    "\n",
    "The derivative of the sigmoid transfer function is:\n",
    "\\begin{align} \n",
    "\\frac{dF}{dx} & = \\frac{d}{dx} (1+\\exp\\{-a(x-\\theta)\\})^{-1}  \\\\\n",
    "& = a\\exp\\{-a(x-\\theta)\\} (1+\\exp\\{-a(x-\\theta)\\})^{-2}. \\qquad (5)\n",
    "\\end{align}\n",
    "\n",
    "Let's now find the expression for the derivative $\\displaystyle{\\frac{dF}{dx}}$ in the following cell and plot it."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "execution": {
     "iopub.execute_input": "2021-05-25T01:14:12.884891Z",
     "iopub.status.busy": "2021-05-25T01:14:12.883728Z",
     "iopub.status.idle": "2021-05-25T01:14:12.885526Z",
     "shell.execute_reply": "2021-05-25T01:14:12.885981Z"
    }
   },
   "outputs": [],
   "source": [
    "def dF(x, a, theta):\n",
    "  \"\"\"\n",
    "  Population activation function.\n",
    "\n",
    "  Args:\n",
    "    x     : the population input\n",
    "    a     : the gain of the function\n",
    "    theta : the threshold of the function\n",
    "\n",
    "  Returns:\n",
    "    dFdx  : the population activation response F(x) for input x\n",
    "  \"\"\"\n",
    "\n",
    "  ###########################################################################\n",
    "  # TODO for students: compute dFdx ##\n",
    "  raise NotImplementedError(\"Student excercise: compute the deravitive of F\")\n",
    "  ###########################################################################\n",
    "\n",
    "  # Calculate the population activation\n",
    "  dFdx = ...\n",
    "\n",
    "  return dFdx\n",
    "\n",
    "\n",
    "pars = default_pars_single()  # get default parameters\n",
    "x = np.arange(0, 10, .1)      # set the range of input\n",
    "\n",
    "# Uncomment below to test your function\n",
    "# df = dF(x, pars['a'], pars['theta'])\n",
    "# plot_dFdt(x, df)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 433
    },
    "colab_type": "text",
    "execution": {
     "iopub.execute_input": "2021-05-25T01:14:12.906538Z",
     "iopub.status.busy": "2021-05-25T01:14:12.904302Z",
     "iopub.status.idle": "2021-05-25T01:14:13.113944Z",
     "shell.execute_reply": "2021-05-25T01:14:13.113419Z"
    },
    "outputId": "540be28c-ffc2-44d3-fbdb-6f27b3305405"
   },
   "source": [
    "[*Click for solution*](https://github.com/NeuromatchAcademy/course-content/tree/master//tutorials/W2D4_DynamicNetworks/solutions/W2D4_Tutorial1_Solution_ce2e3bc5.py)\n",
    "\n",
    "*Example output:*\n",
    "\n",
    "<img alt='Solution hint' align='left' width=560 height=416 src=https://raw.githubusercontent.com/NeuromatchAcademy/course-content/master/tutorials/W2D4_DynamicNetworks/static/W2D4_Tutorial1_Solution_ce2e3bc5_0.png>\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text"
   },
   "source": [
    "## Exercise 5: Compute eigenvalues\n",
    "\n",
    "As discussed above, for the case with $w=5.0$ and $I_{\\text{ext}}=0.5$, the system displays **three** fixed points. However, when we simulated the dynamics and varied the initial conditions $r_{\\rm init}$, we could only obtain **two** steady states. In this exercise, we will now check the stability of each of the three fixed points by calculating the corresponding eigenvalues with the function `eig_single`. Check the sign of each eigenvalue (i.e., stability of each fixed point). How many of the fixed points are stable?\n",
    "\n",
    "Note that the expression of the eigenvalue at fixed point $r^*$\n",
    "$$\\lambda = [-1+ wF'(w\\cdot r^* + I_{\\text{ext}};a,\\theta)]/\\tau$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "execution": {
     "iopub.execute_input": "2021-05-25T01:14:13.121307Z",
     "iopub.status.busy": "2021-05-25T01:14:13.120778Z",
     "iopub.status.idle": "2021-05-25T01:14:13.124337Z",
     "shell.execute_reply": "2021-05-25T01:14:13.123820Z"
    }
   },
   "outputs": [],
   "source": [
    "def eig_single(fp, tau, a, theta, w, I_ext, **other_pars):\n",
    "  \"\"\"\n",
    "  Args:\n",
    "    fp   : fixed point r_fp\n",
    "    tau, a, theta, w, I_ext : Simulation parameters\n",
    "\n",
    "  Returns:\n",
    "    eig : eigevalue of the linearized system\n",
    "  \"\"\"\n",
    "  #####################################################################\n",
    "  ## TODO for students: compute eigenvalue and disable the error\n",
    "  raise NotImplementedError(\"Student excercise: compute the eigenvalue\")\n",
    "  ######################################################################\n",
    "  # Compute the eigenvalue\n",
    "  eig = ...\n",
    "\n",
    "  return eig\n",
    "\n",
    "\n",
    "# Find the eigenvalues for all fixed points of Exercise 2\n",
    "pars = default_pars_single(w=5, I_ext=.5)\n",
    "r_guess_vector = [0, .4, .9]\n",
    "x_fp = my_fp_finder(pars, r_guess_vector)\n",
    "\n",
    "# Uncomment below lines after completing the eig_single function.\n",
    "\n",
    "# for fp in x_fp:\n",
    "#   eig_fp = eig_single(fp, **pars)\n",
    "#   print(f'Fixed point1 at {fp:.3f} with Eigenvalue={eig_fp:.3f}')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text"
   },
   "source": [
    "**SAMPLE OUTPUT**\n",
    "\n",
    "```\n",
    "Fixed point1 at 0.042 with Eigenvalue=-0.583\n",
    "Fixed point2 at 0.447 with Eigenvalue=0.498\n",
    "Fixed point3 at 0.900 with Eigenvalue=-0.626\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 70
    },
    "colab_type": "text",
    "execution": {
     "iopub.execute_input": "2021-05-25T01:14:13.130735Z",
     "iopub.status.busy": "2021-05-25T01:14:13.130223Z",
     "iopub.status.idle": "2021-05-25T01:14:13.134428Z",
     "shell.execute_reply": "2021-05-25T01:14:13.133985Z"
    },
    "outputId": "70e87372-6a77-4376-93d5-a2de111ec1ad"
   },
   "source": [
    "[*Click for solution*](https://github.com/NeuromatchAcademy/course-content/tree/master//tutorials/W2D4_DynamicNetworks/solutions/W2D4_Tutorial1_Solution_e285f60d.py)\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text"
   },
   "source": [
    "## Think! \n",
    "Throughout the tutorial, we have assumed $w> 0 $, i.e., we considered a single population of **excitatory** neurons. What do you think will be the behavior of a population of inhibitory neurons, i.e., where $w> 0$ is replaced by $w< 0$? "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab": {},
    "colab_type": "text",
    "execution": {
     "iopub.execute_input": "2021-05-25T01:14:13.137994Z",
     "iopub.status.busy": "2021-05-25T01:14:13.137488Z",
     "iopub.status.idle": "2021-05-25T01:14:13.142375Z",
     "shell.execute_reply": "2021-05-25T01:14:13.141932Z"
    }
   },
   "source": [
    "[*Click for solution*](https://github.com/NeuromatchAcademy/course-content/tree/master//tutorials/W2D4_DynamicNetworks/solutions/W2D4_Tutorial1_Solution_579bc9c9.py)\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text"
   },
   "source": [
    "---\n",
    "# Bonus 2: Noisy input drives the transition between two stable states\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text"
   },
   "source": [
    "## Ornstein-Uhlenbeck (OU) process\n",
    "\n",
    "As discussed in several previous tutorials, the OU process is usually used to generate a noisy input into the neuron. The OU input $\\eta(t)$ follows: \n",
    "\n",
    "$$\\tau_\\eta \\frac{d}{dt}\\eta(t) = -\\eta (t) + \\sigma_\\eta\\sqrt{2\\tau_\\eta}\\xi(t)$$\n",
    "\n",
    "Execute the following function `my_OU(pars, sig, myseed=False)` to generate an OU process."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "cellView": "form",
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 286
    },
    "colab_type": "code",
    "execution": {
     "iopub.execute_input": "2021-05-25T01:14:13.169517Z",
     "iopub.status.busy": "2021-05-25T01:14:13.164614Z",
     "iopub.status.idle": "2021-05-25T01:14:13.346803Z",
     "shell.execute_reply": "2021-05-25T01:14:13.345925Z"
    },
    "outputId": "5fcb932d-ba75-4b9f-8bb4-030363a53dd7"
   },
   "outputs": [],
   "source": [
    "# @title OU process `my_OU(pars, sig, myseed=False)`\n",
    "\n",
    "# @markdown Make sure you execute this cell to visualize the noise!\n",
    "\n",
    "\n",
    "def my_OU(pars, sig, myseed=False):\n",
    "  \"\"\"\n",
    "  A functions that generates Ornstein-Uhlenback process\n",
    "\n",
    "  Args:\n",
    "    pars       : parameter dictionary\n",
    "    sig        : noise amplitute\n",
    "    myseed     : random seed. int or boolean\n",
    "\n",
    "  Returns:\n",
    "    I          : Ornstein-Uhlenbeck input current\n",
    "  \"\"\"\n",
    "\n",
    "  # Retrieve simulation parameters\n",
    "  dt, range_t = pars['dt'], pars['range_t']\n",
    "  Lt = range_t.size\n",
    "  tau_ou = pars['tau_ou']  # [ms]\n",
    "\n",
    "  # set random seed\n",
    "  if myseed:\n",
    "      np.random.seed(seed=myseed)\n",
    "  else:\n",
    "      np.random.seed()\n",
    "\n",
    "  # Initialize\n",
    "  noise = np.random.randn(Lt)\n",
    "  I_ou = np.zeros(Lt)\n",
    "  I_ou[0] = noise[0] * sig\n",
    "\n",
    "  # generate OU\n",
    "  for it in range(Lt - 1):\n",
    "    I_ou[it + 1] = (I_ou[it]\n",
    "                    + dt / tau_ou * (0. - I_ou[it])\n",
    "                    + np.sqrt(2 * dt / tau_ou) * sig * noise[it + 1])\n",
    "\n",
    "  return I_ou\n",
    "\n",
    "\n",
    "pars = default_pars_single(T=100)\n",
    "pars['tau_ou'] = 1.  # [ms]\n",
    "sig_ou = 0.1\n",
    "I_ou = my_OU(pars, sig=sig_ou, myseed=2020)\n",
    "plt.figure(figsize=(10, 4))\n",
    "plt.plot(pars['range_t'], I_ou, 'r')\n",
    "plt.xlabel('t (ms)')\n",
    "plt.ylabel(r'$I_{\\mathrm{OU}}$')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text"
   },
   "source": [
    "## Example: Up-Down transition\n",
    "\n",
    "In the presence of two or more fixed points, noisy inputs can drive a transition between the fixed points! Here, we stimulate an E population for 1,000 ms applying OU inputs."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "cellView": "form",
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 286
    },
    "colab_type": "code",
    "execution": {
     "iopub.execute_input": "2021-05-25T01:14:13.352151Z",
     "iopub.status.busy": "2021-05-25T01:14:13.351593Z",
     "iopub.status.idle": "2021-05-25T01:14:13.609926Z",
     "shell.execute_reply": "2021-05-25T01:14:13.610389Z"
    },
    "outputId": "e327a1fd-2d8a-41dc-a2d9-aef33a54d220"
   },
   "outputs": [],
   "source": [
    "# @title Simulation of an E population with OU inputs\n",
    "\n",
    "# @markdown Make sure you execute this cell to spot the Up-Down states!\n",
    "\n",
    "pars = default_pars_single(T=1000)\n",
    "pars['w'] = 5.0\n",
    "sig_ou = 0.7\n",
    "pars['tau_ou'] = 1.  # [ms]\n",
    "pars['I_ext'] = 0.56 + my_OU(pars, sig=sig_ou, myseed=2020)\n",
    "\n",
    "r = simulate_single(pars)\n",
    "\n",
    "plt.figure(figsize=(10, 4))\n",
    "plt.plot(pars['range_t'], r, 'b', alpha=0.8)\n",
    "plt.xlabel('t (ms)')\n",
    "plt.ylabel(r'$r(t)$')\n",
    "plt.show()"
   ]
  }
 ],
 "metadata": {
  "colab": {
   "collapsed_sections": [],
   "include_colab_link": true,
   "name": "W2D4_Tutorial1",
   "provenance": [],
   "toc_visible": true
  },
  "kernel": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "kernelspec": {
   "display_name": "Python 3",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}
